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AN  OPERATIONAL  UPPER  AIR  ANALYSIS 
USING  THE  VARIATIONAL  METHOD 


Technical  Note  No.  72-3. 
May  1972 


by 


J.  M.  Lewis 


FOREWORD 


This  technical  note  describes  the  Numerical  Variational 
Analysis  (NVA)  program  used  at  FNWC  Monterey  to  prepare  upper 
air  analyses  of  wind  and  temperature  for  the  Global  Band  Grid. 
The  grid  has  a  mesh  length  of  2.5  degrees  at  the  Equator  and 
extends  around  the  world  covering  latitudes  between  40S  and 
60N.  (The  tropical  grid  used  at  FWC  Pearl  Harbor  is  a  sub¬ 
set  of  this  one.) 

The  program  is  unique  to  Naval  Weather  Service  analysis 
effort  in  the  tropics  in  that  it  couples  wind  and  temperature 
through  thermal  wind  constraints  thereby  automatically  ensuring 
vertical  consistency.  It  is  designed  to  make  maximum  use  of 
all  source  data- -especially  over  sparsely  sampled  ocean  areas. 

Reviewed  and  approved  on  20  May  1972. 


W.  S.  HOUSTON,  Jr. 
Captain,  U.  S.  Navy 
Commanding  Officer 


Fleet  Numerical  Weather  Central 
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ABSTRACT 


An  upper  air  objective  analysis  of  wind  and  temperature 
has  been  developed  for  operational  use  at  Fleet  Numerical 
Weather  Central.  The  primary  feature  of  this  analysis  scheme 
is  the  vertical  and  horizontal  extrapolation  of  information 
into  the  data-sparse  regions  from  the  data-rich  surface  and 
jet  aircraft  level.  The  extrapolation  is  based  on  Sasaki’s 
recent  extensions  of  the  variational  analysis  method.  The 
hydrostatic  equation  and  horizontal  momentum  equations  are 
used  as  dynamical  constraints  and  these  equations  are  combined 
to  relate  temperature  and  wind  through  a  generalized  thermal 
wind  relationship.  The  governing  analysis  equations  are 
elliptic  and  amenable  to  solution  by  the  relaxation  process. 
Vertical  coupling  of  both  wind  and  temperature  is  demonstrated 
through  the  detailed  examination  of  a  case  study  in  the  North 
Pacific. 


vi 


1.  INTRODUCTION 

The  synoptic  analysis  of  meteorological  fields  at 
mid- tropospheric  levels  has  been  plagued  by  the  scarcity 
of  conventional  data.  The  data . coverage  over  the  densely 
populated  continents  of  Europe  and  North  America  has  general¬ 
ly  been  adequate  for  the  definition  of  the  large-scale  or 
extra-tropical  cyclone  features.  For  ocean  operations, 
however,  the  conventional  data  is  strikingly  inadequate  to  . 
define  even  these  large-scale  features.  To  augment  the  . 
information  at  these  mid-levels,  the  relatively  rich  sources 
of  data  at  the  surface  have  been  used..  Yanai  (1964)  achieved 
a.  considerable  improvement  over  standard  methods  when  he 
relied  heavily  upon  the  surface  data,. to  build  upper-air 
analyses  of  wind  in  the  Caribbean.  .  His  success  stemmed  from 
the  ability  to  construct  an  accurate  surface  wind  analysis 
from  relatively  abundant  ship  and  island  reports.  Adjoining 
upper  air  charts  were  then  built  stepwise  by  differential 
analysis  using  the  successive  corrections  method  (SCM) ,  an 
operationally  attractive  scheme  formulated  by  Bergthbrsson 
and  Doos  (1955)  and  modified  by  Cressman  (1959). 

Another  approach  that  has  been  successful  at  the  National 
Hurricane  Research  Laboratory  (NHRL)  is  the  layer  mean  analy¬ 
sis.  In  this  analysis,  the  data  base  for  the  lower  tropo¬ 
sphere  (1000  -  600  mb)  and  the  upper  troposphere  (600  -  200  mb) 
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is  effectively  increased  by  amalgamating  the  information 
from  the  individual  mandatory  levels.  The  apparent  success 
of  this  scheme  is  based  on  the  two-layer  structure  of  the 
tropics  during  most  of  the  year  (Wise  and  Simpson,  1971). 

In  addition  to  increasing  the  data  base  for  each 
meteorological1 variable ,  the  ability  to  couple  variables 
through  dynamical  or  empirical  relations  is  germane  to 
analysis  in  data-sparse  regions.  Thus,  the  observation  of 
a  particular  variable  can  indirectly  produce  information 
about  another  variable.  Sasaki  (1958,  1969)  has  developed 
a  scheme  based  on  variational  calculus  which  couples  the 
meteorological  fields  by  constraining  the  analysis  to 
satisfy  a  governing  set  of  diagnostic  or  prognostic  equations. 
This  methodology  has  been  used  at  Fleet  Numerical  Weather 
Central  (FNWC) "to  produce  surface  winds  and  pressure  on  a 
glbbal-band  grid  which  encompasses  the  tropical  belt  (Lewis 
and  Grayson,  1972) .  The  work  of  Yanai  (1964)  also  coupled 
the  fields  at  the  surface  by  using  a  statist ical -dynamical 
model  proposed  for  use  at  the  Japan  Meteorological  Agency 
(JMA)  and  reported  by  Masuda  and  Arakawa  (1962). 

The  coupling  of  data  in  upper-air  analysis  has  received 
considerable  attention  in  recent  years.  This  interest  has 
been  stimulated  by  the  necessity  to  define  observational 
networks  for  numerical  weather  prediction  (Doos,  1970).  The 
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operational  centers  have  traditionally  used  geostrophic 
theory  in  conjunction  with  the  successive  corrections  method. 
V,rith  upsurge  in  global  forecasting,  however,  the  analysis 
techniques  have  incorporated  ageostrophic  effects.  The 
National  Meteorological  Center  (NMC)  no w  analyzes  heights 
and  winds  by  spectral  methods  which  use  Hough  functions 
(eigenfunctions  of  Laplace's  tidal  equation)  for  horizontal 
modes  and  empirical  orthogonal  functions  for  the  vertical 
modes  (Flattery,  1971)  .  The  results  indicate  that  reasonable 
height  and  wind  fields  are  obtained,  even  at  the  equator. 

Since  the  Hough  functions  are  derived  from  equations  which 
neglect  the  nonlinear  advection  terms  in  the  momentum  equations, 
there  have  been  some  problems  when  the  analysis  is  used  for 
initializing  the.  primitive  equation  forecast  model  (Flattery, 
1970).  Gradient  wind  corrections  have  been  made  in  an 
attempt  to  account  for  the  nonlinear  imbalance. 

The  proposed  upper -air  analysis  scheme  at  FNWC  is  designed 
for  operational  use  on  a  global  band  grid  extending  from 
40S  to  60N  which  has  a  mesh  length  of  150  n  mi  at  the  equator 
and  75  n  mi  at  60N.  The  change  of  mesh  size  is  a  consequence 
of  the  distortion  inherent  in  the  Mercator  secant  projection 
(true  at  22. 5N  and  22. 5 S )  used  for  the  analysis.  The  flight 
support  for  turbo-prop  aircraft  has  been  a  primary  considera¬ 
tion  in  the  development  of  this  scheme.  These  planes  cruise 
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in  the  vicinity  of  the  400  mb  surface  and  flight  support 
and  guidance  at  this  level  are  necessary.  As  mentioned 
previously,  the  data-rich  levels  are  the  surface  and  the  jet 
aircraft  level  (vicinity  of  250  mb).  The  extrapolation  of 
data  from  these  data-abundant  regions  into  the  mid- troposphere 
is  accomplished  by  using  Sasaki's  numerical  variational 
analysis  (NVA)  method.  This  approach  reduces  the  analysis 
problem  to  a  boundary  value  problem  where  the  upper  and  lower 
surface  are  specified  and  the  interior  is  determined  subject 
to  the  governing  elliptic  analysis  equations.  Since  there 
are  scattered  observations  between  the  surface  and  250  mb, 
the  analysis  scheme  must  have  the  capability  of  weighting 
this  local  information  as  well  as  incorporating  the  extrapolated 
information.  The  variational  analysis  scheme  provides  a 
natural  way  to  accommodate  both  the  observations  and  the 
dynamically  extrapolated  information. 

The  surface  analysis  which  is  used  to  provide  boundary 
values  for  the  upper-air  analysis  is  based  on  the  operational 
scheme  reported  by  Lewis  and  Grayson  (1972).  This  analysis 
has  been  subjected  to  a  two-dimensional  adjustment  by  NVA. 

The  momentum  equations  used  as  constraints  for  this  surface 
analysis  are  consistent  with  the  upper-air  scheme  except 
that  they  contain  the  effects  of  surface  friction.  The 
250  mb  analysis  is  found  by  the  SCM.  A  two-dimensional 


adjustment  could  be  performed  at  this  level,  but  has  been 
postponed  until  the  general  features  of  the  extrapolation 
procedure  are  clarified.  In  order  to  incorporate  dynamical 
information  at  this  level,  the  forecasted  wind  from  the  FNWC 
primitive  equation  model  (Kesel  and  Winninghoff,  1971)  is 
used  as  the  guess  field  for  the  SCM  between  20N  -  60N. 

Since  aircraft  rarely  report  geopotential  (D-values)  and, 
moreover,  because  of  the  limited  credibility  of  these  observa¬ 
tions,'  the  coupling  of  temperature  and  wind  has  been  used 
to  build  the  upper-air  analyses.  The  subjective  or  manual 
analysis  of  upper-air  wind  has  heavily  relied  upon  this 
coupling  or  thermal  wind  relationship  (see  Saucier,  1955). 

In  order  to  couple  these  fields  over  the  tropics  and  extra- 
tropical  regions,  a  more  general  form  of  the  thermal  wind 
relation  which  includes  ageostrophic  effects  must  be  employed. 
Forsythe  (1945)  developed  graphical  techniques  which  enabled 
the  hand-analyst  to  account  for  ageostrophic  effects  between 
the  temperature  gradient  and  wind  shear.  It  is  precisely 
this  traditional  synoptic  rule  connecting  temperature  gradient 
and  wind  shear  that  provides  the  vehicle  for  extrapolation 
in  the  objective  scheme. 
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2.  BASIC  INGREDIENTS  FOR  THIS  ANALYSIS,  OBSERVATIONS 

AND  THE  GENERAL  THERMAL  WIND  RELATION 

The  operational  approach  to  numerical  variational 
analysis  (referred  to  as  NVA  in  future  discussion)  adopted 
at  FNWC  is:  (1)  irregularly  distributed  observations  of 
each  variable  are  independently  interpolated  to  grid  points 
by  the  SCM,  (2)  a  reliability  is  attached  to  the  variable 
at  each  grid  point  and  is  dependent  upon  both  the  number 
of  observations  affecting  this  point  and  the  distance  of  the 
observations  from  the  point,  and  (3)  an  adjustment  of  the 
"observations"  (the  analysis  resulting  from  the  SCM)  by  the 
NVA  equations  which  have  incorporated  the  model  constraints. 

The  analysis  is  performed  every  6  hr  and  the  natural 
candidate  for  the  guess  field  of  SCM  is  the  previous  analysis 
or,  equivalently  stated,  the  6  hr  persistence  field.  In 
addition  to  the  persistence  field,  however,  the  FNWC  primitive 
equation  model  produces  a  forecast  every  12  hr  over  the 
Northern  Hemisphere.  The  guess  field  for  the  wind  is  6  hr 
persistence  south  of  20N.  Between  20N  -  60N,  the  6  hr  and 
12  hr  forecasts  are  alternately  used.  That  is,  the  06  GMT 
and  18  GMT  global-band  analyses  use  the  6  hr  forecasts  from 
00  GMT  and  12  GMT,  respectively.  Similarly,  the  00  GMT 
and  12  GMT  analyses  use  12  hr  forecasts  from  12  GMT  (previous 
day)  and  00  GMT,  respectively.  The  temperature  guess  field, 
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however,  is  strictly  the  6  hr  persistence  field.  The  feed¬ 
back  from  wind  to  temperature  occurs  in  the  NVA  scheme. 

Only  when  the  "observations"  (SCM  analyses)  are  independent 
can  the  variational  adjustment  process  guarantee  maximum 
likelihood  estimates  (Cramer,  1955;  Lewis  and  Grayson,  1972). 

Four  successive  scans  are  made  to  correct  the  guess 
field.  For  the  analysis  of  temperature,  the  radii  are  600, 

300,  200  and  100  n  mi,  respectively.  Grid  points  which  are 
600  n  mi  or  further  from  a  report  are  unmodified,  i.e., 
the  6  hr  persistence  becomes  the  current  analysis.  It  is 
precisely  this  kind  of  event  that  often  causes  the  SCM 
to  produce  "unmeteorological"  fields,  i.e.,  the  data-sparse 
areas  adjoining  the  data-rich  areas  show  discontinuities 
which  are  characteristic  of  the  mathematical  weighting  func¬ 
tions  of  the  SCM  and  are  not  attributable  to  meteorological 
events. 

In  a  similar  fashion,  three  scans  are  used  by  the  SCM 
for  wind  analysis.  Here  the  two  component  directions  are 
analyzed  together.  In  addition  to  speeding  the  analysis 
procedure,  the  decision  to  accept  or  reject  wind  data  is 
based  on  the  vector  rather  than  the  separate  scaler  components. 
The  radii  for  the  wind  analysis  are  450,  300  and  150  n  mi, 
respectively.  The  successively  decreasing  scan  radii  are 
designed  to  capture  the  larger  scale  features  on  the  preliminary 
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scan  and  then  include  the  smaller  scale  features  on  the 
following  scans.  Stephens  (1967)  has  discussed  the  relation¬ 
ship  between  characteristic  meteorological  scale  and  the  scan 
radii  under  the  assumption  of  uniformly  distributed  data. 

The  second  ingredient  to  NVA  is  the  model  dynamics 
which  will  be  used  to  couple  the  variables  and  provide  the 
path  to  extrapolation  from  the  data-dense  to  data-sparse 
regions.  The  governing  equations  are  modified  in  a  manner 
consistent  with  the  NVA  at  the  surface  (Lewis  and  Grayson, 
1972).  Namely,  the  nonlinear  advection  terms  are  approximated 
by  the  distributions  resulting  from  the  SCM.  The  primary 
advantage  is  the  simplification  of  the  Euler  equations  and 
the  guarantee  of  an  elliptic  differential  equation  system 
regardless  of  the  weights  used  in  the  variational  formulation. 
The  guaranteed  ellipticity  leads  to  rapid  convergence  when 
the  relaxation  process  is  used.  The  governing  equations  in 
Mercator  coordinates  are: 


3 1 

V*  Vu  - 

fv  =  -m 

3  <p 

3x 

(1) 

+ 

Jt 

a  •  w  + 

fu  =  -m 

3<J> 

w 

(2) 

= 

3a 

RT 

(3) 

a  =  In 

(P0'P> 

(4) 
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where:  t 


time 


east-west  coordinate  (positive  eastward) 
north-south  coordinate  (positive  northward) 
vertical  coordinate  (positive  upwards) 


=  ui  +  vj 


Ul  +  V] 


horizontal  velocity- 
unit  vector  along  positive  x 
unit  vector  along  positive  y 
velocity  from  SCM 


temperature 
geopotential 
Coriolis  parameter 

cos0Q  sec0  :  image  scale  factor  for  Mercator 

projection  true  at  22. 5N  and  22.5S 


latitude  (0  =  22.5  degrees) 
reference  pressure  (1013  mb) 
gas  constant  for  dry  air 

m  i  +  m  ~  j  :  gradient 


9y 


operator 


To  simplify  the  notation,  the  nonlinear  advection  terms 


are  denoted  by: 


A  =  \V»  Vu 


B  =  \V  •  Vv 


The  general  thermal  wind  relation  is  derived  by  eliminating 
the  geopotential  from  the  momentum  equations,  (1)  and  (2), 
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by  substituting  from  the  hydrostatic  equation,  (3).  In 
component  form,  this  relation  is 


3 

Jt 

(— i  = 
'■3aJ 

3 

'  'So 

A  +  f 

3v 

3a 

mR 

3T  = 
3x 

3 

3t 

f— i  = 

3 

3a 

B  -  f 

3u 

So  " 

mR 

3T 

^y  = 

(7) 

(8) 


where  the  separation  of  terms  by  functions  F  and  G  will  prove 
useful  in  the  discussion  on  minimization. 
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3.  THE  VARIATIONAL  FORMULATION  AND  ASSOCIATED  EULER  EQUATIONS 

The  timewise  localized  version  of  NVA  (Sasaki,  1970)  is 
adopted  for  the  upper-air  analysis  and  the  functional  to  be 
minimized  assumes  the  form: 


where 

a  :  Gauss'  precision  modulus  for  wind 

1  . 

3  :  Gauss'  precision  modulus  for  temperature 

a  :  Dynamical  weight 

x  :  east-west  coordinate  on  earth 
e 

y  :  north-south  coordinate  on  earth 
(~)  :  analysis  resulting  from  the  SCM 

xe'ye 

J( - — )  :  Jacobian  of  transformation  from  earth  to 

*  * » y 

Mercator  coordinates 
dV  :  incremental  volume  (dxdydo) 

The  relationship  between  earth  and  map  coordinates  is 

X  fy 

dx£dye  =  dxdy  •  (10) 

Consequently,  (9)  implies  that  equal  weight  is  given  to 
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equal  earth  areas,  not  equal  map  areas.  For  the  Mercator 
pro j  ection , 

■  '  x  ,y  „ 

..  J(-hr)  ■ m  .(11) 

By  virtue  of  (7)  and  (8),  the  functional  can  be  rewritten 


I  = 


r  -  2  2  ^ 
a[(u-u)z  +  (v-vj  ] 

+  3 [ (T-T) 2] 


+  a [F2  +  G2] 


J(- 


x  ,y 
e  * 7  e 

x,y 


)  dV 


(12) 


The  first  set  of  bracketed  terms  force  the  analysis  toward 
the  SCM  velocity  in  proportion  to  a.  Similarly,  the  second 
bracketed  term  forces  the  temperature  toward  the  SCM 
temperature  in  proportion  to  8.  The  weights  5  and  3  are 
called  observational  weights  or  Gauss'  precision  moduli 
and  specified  beforehand,  i.e.,  before  the  minimization 
process.  They  are  a  measure  of  the  confidence  in  the 
respective  SCM  fields.  Following  the  discussion  of  these 
weights  in  the  paper  by  Lewis  and  Grayson  (1972),  let 


5  =  l/crw  ,  5  =  l/oT  (13) 

where  aw  is  the  mean-square  difference  between  the  separate 
wind  components  of  NVA  and  SCM  and  a^,  is  the  mean-square 
difference  between  the  NVA  and  SCM  temperatures. 

In  analogy  with  the  observational  weights  5  and  3,  the 
dynamical  weight  a  is  expressed  as 
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,(14) 


a  =  l/at 

where  is  the  mean-square  value  of  local  time  rate  of 
change  of  wind  shear  as  calculated  from  (7)  and  (8).  Alter¬ 
natively  stated,  a  is  a  measure  of  the  steadiness  or  degree 
with  which  u,  v,  and  T  satisfy  F  =  G  =  0.  As  ct  increases, 
the  analyzed  fields  are  forced  closer  to  the  steady  state 
relationships. 

The  observational  weights  5  and  8  should  be  a  function 
of  space  (and  time)  and  generally  distinguish  between  (1)  grid 
points  that  benefited  from  nearby  observations  and  (2)  grid 
points  unaffected  by  observations  and  consequently  reflecting 
6  hr  persistence  after  the  SCM  analysis.  A  method  for  speci¬ 
fication  of  a  and  8  in  terms  of  the  distance-dependent  weighting 
implicit  in  the  SCM  has  been  discussed  by  Lewis  and  Grayson 
(1972).  During  the  development  stages  of  the  upper  air  analy¬ 
sis,  the  weights  have  not  varied  as  a  function  of  space  but 
have  been  chosen  to  reflect  only  typical  variance  of  the 
wind  and  temperature. 

As  mentioned  in  the  Lewis  and  Grayson  (1972)  paper,  the 
separately  specified  values  of  aw,  and  a t  are  not  preserved 
in  the  minimization  process.  However,  the  ratios  between 
these  values  are  preserved.  This  fact  will  be  obvious  when 
the  Euler  equations  (Eqs.  17-19)  are  examined.  If  a,  b  and 
c  are  chosen  for  ow,  and  at,  the  final  analysis  produces 
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values  d,  e  and  f,  respectively,  where  a:b:c  =  d:e:f.  Rather 
than  specifying  a  ,  aT  and  c^,  it  is  sufficient  to  specify 
aw/at  and  o^/o^  *  . 

A  necessary  and  sufficient  condition  for  minimization 
of  (12)  is  the  vanishing  of  the  first  variation,  denoted  by 
61.  The  sufficiency  of  this  condition  is  directly  attributa¬ 
ble  to  the  quadratic  nature  of  the  integrand  (Lanczos,  1970, 
p.  41).  Using  integration  by  parts  and  the  commutative 
properties  of  the  6-process  (Lanczos,  1970,  p.  56),  the 
functional  takes  the  form 

0  ■  61  (p6u  +  +  dV  *  {surfaceVbiundaries}  •  (15) 

'V 

The  expressions  P,  Q  and  W  are  differential  expressions  which 
will  be  explicitly  written  later. 

The  boundary  terms  indicated  on  the  right-hand  side  of 
(15)  involve  integrals  over  the  bounding  surfaces,  viz . , 
the  walls  at  40S  and  60N,  the  250  mb  level  and  the  earth’s 
surface.  The  east-west  boundaries  are  cyclic  because  of  the 
360°  longitudinal  span  of  the  analysis  region.  The  conditions 
necessary  to  cause  these  integrals  to  vanish  are  called  the 
natural  boundary  conditions.  There  is  generally  more  than 
one  set  of  acceptable  boundary  conditions.  The  simplest 
boundary  conditions  are:  6T  =  0  on  the  north  and  south 
walls  and  6u  =  6v  =  0  along  the  upper  and  lower  surfaces. 
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These  are  Dirichlet  conditions  and  allow  no  adjustment  of 
temperature  and  velocity  along  the  boundaries  where  the 
respective  variations  vanish.  Alternative  boundary  condi¬ 
tions  can  be  found  by  examination  of  the  integrals  given 
in  the  appendix.  The  simplicity  of  the  Dirichlet  conditions, 
especially  in  regard  to  the  programming  of  numerical  methods 
for  solution,  has  led  to  the  adoption  of  this  set  of  boundary 
conditions. 

Having  disposed  of  the  terms  evaluated  on  the  surface 
boundaries  and  because  of  the  independence  of  the  variations 
of  u,  v  and  T,  the  solution  to  (15)  is 


p  =  q  =  W  =  0 


(16) 


where 


P  =  •< 


V. 

r 


Q  =  < 


r2  3  U  /Ct\  /  n 

f  rr  •  y  (u'u) 

dC 


CT.  9  r  9  T  \  r  9  B 

+  mfR  3o  +  f  TT 

'  dO  j 


>  =  o 


f2  2-t  -  (f)(v-v) 
BO 


a' 


rn  3  r3T>  .  32A 
•  m£R  -  £  “T 


3o‘ 


r  =  o 


(17) 


(18) 


15 


w  = 


V-(m'2VT)  -  (mR)'2  (|)(T-T) 


ll-V*  (f_  +  V  •(-£-) 

m  W  J<TJ 

m  R 


(19) 


/«  _  3A  *  3B  ^ 
1  ■  -rr  1  +  ) 


(20) 


/V 

The  local  vertical  unit  vector  is  k  and  the  dot  and  cross- 
products  have  been  incorporated  in  (19).  Equations  (17) 
and  (18)  can  be  combined  to  produce  the  vector  representation: 

®  © 


f2  _  (|}  (  vy  . 

3o 


+  fRkx  VT  -  fkx  4-  (C 
da  9a 


(21) 


The  equations,  (17) -(19),  or  equivalently,  (19)  and  (21),  are 
called  The  Euler  equations.  Together  with  the  Dirichlet 
boundary  conditions,  they  constitute  the  governing  equations 
for  the  upper-air  analysis. 

The  derivation  of  (17)- (19)  has  assumed  that  a  is  invariant 
in  space.  A  variable  a  can  easily  be  included  in  the  equations, 
but  the  applicability  of  the  model  constraints  over  the  entire 
analysis  region  dictates  the  constancy.  The  ratios  8/a  and 
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a/a  are  the  controlling  factors  in  (19)  and  (21),  respectively. 
As  the  reliabilities  of  the  SCM  wind  and  temperature  increase, 
the  ratios  a/a  and  3/a  increase,  respectively.  Consequently, 
term  (2)  in  both  (19)  and  (21)  forces  the  analysis  toward  the 
SCM  analysis  in  proportion  to  these  ratios. 

Terms  (T)  and  @  in  (21)  represent  the  geostrophic 
coupling  between  wind  and  temperature.  In  fact,  except  for 
the  differentiation  with  respect  to  a,  this  is  the  geostrophic 
thermal-wind  relationship.  Interpreted  physically,  the  verti¬ 
cal  change  of  horizontal  temperature  gradient  contributes 


2 

r\  *■*  iir 

to  the  curvature,  — ,  of  the  wind  profile.  Term  @ 

9cr 

brings  the  nonlinear  advection  of  momentum  into  the  process 
and  is  generally  an  important  contributer  in  the  regions  of 
strong  horizontal  shear  such  as  in  the  area  neighboring  the 
jet  stream.  Although  this  term  is  approximated  by  the  SCM 
analysis,  the  approximation  is  improved  by  a  "re-cycling" 
technique  explained  in  the  next  chapter. 

2 

The  horizontal  curvature  of  temperature,  V  T,  and  the 
vertical  variation  of  vorticity,  (1c •  V x  \V)  ,  are  related  by 
terms  (T)  and  (3)  in  (19).  Again,  these  terms  stem  from  a 
strict  hydrostatic-geostrophic  balance.  The  nonlinear 
advective  effects  are  contained  in  (3)  and  are  associated 
with  vertical  variation  of  streamline  curvature  as  expressed 

by  U  ^  ( — j) ,  etc.  These  nonlinear  terms  are  subject  to 
9x 
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improved  approximation  by  the  re-cycling  method  mentioned 
above. 

In  the  development  of  the  finite-difference  analogs  of 

the  Euler  equations,  a  departure  from  the  exact  form  has 

been  made.  As  shown  by  Forsythe  and  Wasow  (1960,  p.  182), 

the  construction  of  symmetric  operators  follows  directly 

from  the  finite-difference  form  of  the  functional,  Eq.  (12) 

in  this  case.  Assuming  that  a  centered  difference  is  used 

9  \V 

to  approximate  in  the  functional,  e.g., 


9  W, 

9a  -*400  mb 


(22) 


then  the  finite-difference  analog  to  — should  include 

9a 


values  at  both  the  second  level  above  and  the  second  level 


below  the  one  in  question.  Because  of  the  limited  vertical 
extend  of  the  discretized  analysis  space,  only  the  first 
level  above  and  the  first  level  below  are  used  in  calculating 
the  second  derivative.  This  departure  from  rigorous  form 
introduces  some  discrepancy  between  the  exact  minimum  for 
(12)  and  the  value  obtained  by  using  the  solution  to  the 
assumed  finite-difference  form  of  (17) -(19).  However,  the 


basic  concept  of  coupling  the  variables  and  extrapolation 
remains  intact. 


4.  METHOD  OF  SOLUTION 


The  set  of  Euler  equations  are  coupled  elliptic 
differential  equations  (Helmholtz  type)  subject  to 
Dirichlet  conditions.  The  form  of  the  equation  lends 
itself  to  a  vertically  staggered  grid  arrangement,  viz . , 
an  alternating  sequence  of  temperature  and  wind  fields 
along  the  a-axis.  Equation  (19)  associates  the  vertical 
gradient  of  vorticity,  term  (?)  ,  and  vertical  gradient  of 
momentum  advection,  term  (7)  ,  with  the  horizontal  tempera¬ 
ture  distribution.  Consequently,  the  finite-difference 
g  \V 

evaluation  of  can  be  constructed  using  only  the  levels 
adjacent  to  the  temperature.  A  similar  evaluation  is  possible 
in  connection  with  (21)  and  the  finite-difference  form  of 
A  principal  advantage  of  the  staggering  scheme  is  the 
reduction  in  both  computer  memory  and  time  requirements. 

The  staggered  system  with  the  assignments  for  temperature 
and  wind  is  shown  in  Fig.  1.  The  operational  sequence  of 
events  in  the  analysis  package  is  (1)  the  surface  wind  and 
pressure  field  are  analyzed  by  the  two-dimensional  NVA 
scheme  (Lewis  and  Grayson,  1972),  (2)  the  SCM  upper-air 
analysis  is  used  for  temperature  and  wind  and  (3)  the  varia¬ 
tional  analysis  is  accomplished  through  the  solution  of 
(17)- (19).  The  solution  of  these  Euler  equations  is  schematical 
ly  represented  in  Fig.  1. 
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An  elaboration  upon  the  cycling  process  is  instructive: 


STEP  1. 


STEP  2. 


STEP  3. 


STEP  4. 


The  S50  mb  temperature  is  found  by  solving  (19) 

by  the  accelerated  Liebmann  method  (Miyakoda,  1961). 

On  the  first  cycle,  the  available  wind  field  at  the 

adjacent  upper  level  is  the  SCM  analysis. 

The  700  mb  wind  field  is  found  by  solving  (17)  and 

(18) ,  again  by  the  accelerated  Liebmann  method. 

The  temperature  field  at  850  mb  found  in  STEP  1  is 

now  used  in  conjunction  with  the  SCM  temperature 

3T 

analysis  at  500  mb  to  calculate  .  The  wind  fields 

da 

at  400  mb  (SCM  on  this  first  cycle)  and  the  surface 

2 

3  41  \V 

are  used  to  estimate  — . 

3a 

The  temperature  at  500  mb  is  found  by  solving  (19), 
which  uses  NVA  fields  below  this  level  and  SCM 
fields  above. 

The  w*ind  field  at  400  mb  is  found  by  solving  (17) 
and  (18),  which  use  NVA  fields  below  and  SCM  fields 


above. 

STEP  5.  The  temperature  at  300  mb  is  found  by  solving  (19), 
which  uses  the  NVA  fields  at  all  levels  except  250 
mb  where  the  wind  analysis  is  derived  by  the  SCM 
scheme. 

With  the  completion  of  this  complete  cycle  through  the 
a-layers,  a  re-cycling  process  is  performed.  Now,  the  adjusted 
fields  from  the  first  cycle  are  used  as  forcing  functions  in 
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(17) -(19).  As  expected,  the  convergence  of  the  Liebmann 
method  is  much  faster  during  the  re-cycling  (see  Table  3). 

The  nonlinear  advection  terms,  A  and  B,  are  also  re-calculated 
and.  these  improved  estimations  are  used  in  (17)-(19).  The 
process  of  re-cycling  could  be  continued  indefinitely  but 
obviously  at  the  expense  of  computer  time. 

One  re-cycle  through  the  levels  accounts  for  most  of 
the  adjustment  which  is  fortunate  from  an  operational 
viewpoint.  A  more  rigorous  criterion  for  termination  could 
be  based  upon  the  ratios  of  observational  to  dynamical  weight. 
Thus,  having  specified  ow/at  and  Oj/vt  beforehand,  the  final 
analysis  should  approach  these  ratios  as  the  number  of  cycles 
increases . 
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5.  APPLICATION 


The  analysis  scheme  has  been  tested  using  data  collected 
on  15  February  1972.  Although  the  objective  analysis  was 
performed  over  the  entire  grid  (-7000  points  at  each  level), 
the  synoptic  discussion  will  concentrate  on  the  Pacific  Ocean 
region  between  20  and  60N.  The  jet-aircraft  flights  between 
the  Pacific  islands  and  mainland  (Asia  and  North  America) 
provide  credible  winds  for  the  upper  boundary  condition. 
Additionally,  the  ship  traffic  and  island  density  anchor 
the  surface  analysis. 

Figure  2  depicts  the  sea  level  pressure  analysis  via 
two-dimensional  NVA  (Lewis  and  Grayson,  1972).  The  data 
coverage  is  uniformly  dense  and  both  pressure  and  wind 
reports  have  been  superimposed  over  the  pressure  pattern. 

The  surface  wind  analysis  obtained  simultaneously  by  the 
NVA  scheme  is  shown  in  Fig.  3.  This  analysis  serves  as  the 
lower  boundary  condition  for  the  construction  of  the  upper 
air  analysis.  The  wind  has  been  expressed  in  knots  to 
conform  with  the  standard  operational  format. 

The  upper  boundary  condition  is  the  SCM  analysis  of  wind 
at  250  mb.  This  field  is  displayed  in  Fig.  4  and  wind 
observations  incorporated  into  this  analysis  are  shown  in 
Fig.  5.  The  primitive  equation  forecast  (12  hr)  is  used  as 
the  guess  and  accounts  for  the  definition  of  features  in 
data-void  areas.  For  example,  the  pronounced  ridge  near 
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160-170E  and  extending  from  30-45N  is  virtually  the 
unadulterated  guess  field. 

With  the  specification  of  the  boundary  conditions  and  the 
assignment  of  the  weight  ratios  a/a  and  3/a,  the  Euler 
equations  can  be  solved.  The  solution  is  accomplished  by 
using  the  accelerated  Liebmann  method  in  the  level-by- level 
sequence  discussed  in  Chapter  4.  The  convergence  rate  of 
this  iterative  method  is  strongly  dependent  upon  the  weight 
ratios.  As  3/a  and  a/ a  increase,  the  rate  of  convergence 
increases.  Essentially,  term  (T)  in  both  (19)  and  (21) 
becomes  increasingly  important  and  the  coupling  between  the 
wind  and  temperature  is  weakened.  Operationally,  the  range 
of  values  for  these  ratios  must  be  restricted  to  prevent 
the  usage  of  unreasonable  amounts  of  computer  time.  The 
chosen  ratios,  however,  should  guarantee  a  feedback  or  coupling 
between  the  wind  and  temperature.  The  results  discussed  in 
this  paper  have  been  obtained  with  the  following  ratios: 


a/a  =  at/aw  =  1.5 
3/a  =  crt/a,p  =  250 


(23) 


where  the  implied  space,  time  and  temperature  units. are  n  mi, 
hr  and  °C,  respectively. 

In  order  to  attach  physical  meaning  to  these  weights, 
the  root-mean-square  difference  between  the  SCM  fields  and 
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the  adjusted  fields  must  be  examined.  Representative  values 
at  500  mb  and  700  mb  are  displayed  in  Table  1.  It  is  signifi¬ 
cant  from  an  operational  viewpoint  that  the  major  adjustment 
occurs  within  two  complete  cycles  through  the  a-layers.  This 
table  also  depicts  the  root-mean- square  value  of  the  shear 
tendency,  S(F)  and  S(G).  Since  A a  ~  0.9  for  the  layer  between 
the  surface  and  400  mb,  values  of  S(F)  and  S(G)  can  be 
roughly  interpreted  as  the  400-1000  mb  wind  shear  change  over 
a  one  hr  period. 

A  measure  of  the  overall  convergence  (distinguished  from 
the  convergence  at  a  particular  level  on  a  given  cycle)  can 
be  formulated  in  terms  of  the  difference  between  calculated 
and  hypothetical  ratios.  From  the  form  of  (12)  and  the 
definitions  for  the  weights  given  in  (13) ,  the  following 
relationship  should  hold: 


=  S2(T) :S2(u) :S2(v) :S2(F) :S2(G) 


(24) 


where  S  is 

get 


the  root-mean- square  designator. 


r 

$/a  < 


[S(F)/S(T)]2 
[S(G)/S(T) ] 2 


From  (24)  we 


(25) 
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and 


=  [S(F)/S(u) ] 2 
=  [S  (F)  / S  (v)  ]  2 

S/a  J  -  (26) 

=  -[S(G)/S(u)]^ 

=  (S  (G)  /S  (v)  ]  2 

V. 

Average  values  of  a/a  and  8/a  based  on  (25)  and  (26) 
are  displayed  in  Table  2.  The  marked  improvement  within 
two  cycles  is  evident.  However,  the  theoretical  values  can 
never  be  obtained  because  of  the  truncation  error  introduced 
by  the  finite  difference  analogs  to  (17)- (19).  Additionally, 
the  license  used  to  formulate  the  second  derivative  in  terms 
of  adjacent  levels  (see  end  of  Chapter  3)  ultimately  causes 
some  discrepancy. 

As  expected  a  priori ,  the  amount  of  computer  time  necessary, 
to  solve  for  u,  v  or  T  at  a  particular  level  decreases  with 
every  cycle  through  the  o-layers.  Table  3  shows  the  number 
of  iterations  for  both  temperature  and  wind  at  500  and  700  mb, 
respectively.  The  convergence  criterion  was  established  to 
yield  solutions  accurate  to  10  ^  -  10  ^°C  and  10  ^  -  1  kt 
for  temperature  and  wind,  respectively.  These  values  are 
substantially  below  the  accuracy  of  the  observations. 

In  order  to  reveal  the  magnitude  of  adjustment  in  the 
wind  field,  a  vector  difference  between  the  NVA  and  SCM 
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analyses'  is  displayed.  Figures  6-7  show  the  vector  difference 
at  700  and  400  mb,  respectively,  after  two  complete  cycles 
through  the  a- layers.  The  SCM  analysis  at  700  and  400  mb 
is  shown  in  Figs.  8  and  9,  respectively.  An  encouraging 
result  that  appears  on  both  Figs.  6  and  7  is  the  organization 
of  the  adjustment.  That  is,  the  scale  of  the  adjustment  is 
the  same  order  of  magnitude  as  the  synoptic  scale  disturbances 
and  does  not  exhibit  two  or  three  mesh-length  oscillations. 

The  maximum  adjustment  at  700  mb  occurred  in  connection 
with  the  cyclone  at  approximately  40N,  175W.  The  SCM 
analysis  shows  definite  cyclonic  circulation  above  the  surface 
low  pressure  center,  but  the  adjustment  basically  increases 
the  intensity  of  this  vorticity  center.  Whereas  the  SCM 
analysis  indicates  wind  speeds  the  order  of  10  kt  to  the 
west  of  the  center,  the  NVA  scheme  has  speeds  around  30  kt. 
There  were  no  700  mb  observations  in  this  area  to  verify  the 
adjustment,  but  the  strength  of  the  surface  system  lends 
support  to  both  the  sense  and  magnitude  of  the  modification. 
Another  interesting  adjustment  occurred  east  of  the  Japanese 
Islands  in  the  vicinity  of  the  col  in  the  sea  level  pressure. 
The  magnitude  of  700  mb  wind  above  this  area  is  decreased  by 
NVA  and  the  cyclonic  vorticity  is  weakened.  The  cyclonic 
vorticity  is  increased,  however,  directly  above  the  two 
low  pressure  systems.  Generally,  the  circulations  at  the 
surface  have  been  incorporated  into  the  700  mb  chart. 
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The  scale  of  adjustment  at  400  mb  appears  to  be  larger 
than  at  700  mb  as  expected.  For  example,  the  upper-air 
circulation  accompanying  the  mid-Pacific  low  is  adjusted 
in  the  cyclonic  sense  but  over  a  much  larger  east-west 
extent.  A  conspicuous  amplification  of  the  ridge  along  145W 
can  be  seen  by  examination  of  the  400  mb  adjustment  field. 
Again,  this  seems  to  reflect  the  feedback  of  information 
from  250  mb. 

The  coupling  between  the  temperature  and  wild  fields 
is  displayed  in  Fig.  10.  If  the  model  constraints  were 
strictly  geos trophic-hydrostatic ,  the  direction  of  the 
shear  vector  would  be  along  the  isotherms  and  its  magnitude 
would  be  proportional  to  the  isotherm  gradient.  There  are 
some  noticeably  ageostrophic  regions  especially  in  the 
vicinity  of  weak  thermal  gradients.  The  correlation  between 
the  shear  and  the  thermal  pattern  is  evident  in  the  baro- 
clinic  regions,  e.g.,  the  extensive  SW-NE  flow  from  mid- 
Pacific  to  the  Canadian  border  and  the  SW-NE  flow  just  east 
of  Japan. 

With  the  weight  ratios  given  by  (23),  the  influence 
of  the  wind  field  on  the  temperature  field  is  not  as  con¬ 
spicuous  as  the  reverse  coupling,  i.e.,  the  influence  of 
temperature  on  wind.  This  statement  is  further  substantiated 
by  the  results  in  Table  1  indicating  that  the  NVA  temperature 
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is  very  close,  in  an  rms  sense,  to  the  SCM  temperature. 

Stated  mathematically,  term  (T)  in  (19)  dominates  the  other 
terms  in  this  equation.  In  order  to  isolate  the  influence 
of  wind  on  temperature,  an  experiment  was  conducted  in  which 
3=0.  Thus,  term  (?)  vanishes  in  (19)  and  the  temperature 
at  any  level  is  a  function  of  the  wind  field  and  the  boundary 
values  of  temperature  on  the  north  and  south  walls. 

As  expected,  the  rate  of  convergence  is  slowed  down 
considerably  in  this  case;  in  fact,  the  number  of  iterations 
increased  by  a  factor  of  7.  The  850  mb  temperatures  for 
the  case  3=0  and  the  control  case,  i.e.,  weights  given  by 
(23),  are  shown  in  Figs.  11  and  12,  respectively.  The 
maximum  or  minimum  values  or  "central’'  values  are  in  obvious 
disagreement.  Effectively,  the  3=0  case  has  no  anchoring 
information  on  the  interior.  It  is  very  encouraging,  however, 
to  see  the  coupling  between  isotherm  curvature  and  the  wind 
shear  vector.  Especially  interesting  is  the  cut-off  feature 
in  mid-Pacific. 
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6.  CONCLUDING  REMARKS 

The  primary  results  based  on  this  study  are  (1)  the 
variational  method  can  incorporate  detail  into  the  400-700  mb 
wind  fields  via  vertical  coupling  and  (2)  the  iterative 
technique  for  solving  the  analysis  equations  is  efficient 
enough  to  withstand  :the  rigid  time  requirements  of  operational 
analysis  and  prediction.  Based  on  information  from  the  case 
study,  the  estimated  time  for  solution  to  the  Euler  equations 
is  15  min  [Central  Processor  Unit  (CPU)]  on  the  CDC  6500 
computer.  The  SCM  analysis  takes  approximately  15  min  also, 
yielding  a  combined  time  of  30  min  for  the  entire  analysis. 

The  analysis  will  be  routinely  transmitted  (every  6  hr) 
to  the  outlying  weather  centers  in  June,  1972. 

Considerable  effort  is  currently  concentrated  on  increasing 
the  data-base  for  the  analysis.  In  particular:  (1)  the 
incorporation  of  pibals  (pilot  balloons)  information  will 
ultimately  strengthen  the  low-level  wind  field;  (2)  the 
winds  derived  from  satellite  photographs  will  provide  support 
at  both  the  upper-level  ("250  mb)  and  low-level  ("700  mb); 
and  (3)  temperature  information  from  SIRS  soundings  will  be 
added  to  the  conventional  temperature  reports. 
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TABLE  1.  Root-mean- square  differences  between  adjusted 
fields  (NVA)  and  SCM  fields  denoted  by  S(T),  S(u)  and  S(v) 
where  the  differences  are  calculated  at  500  mb  and  700  mb 
for  the  temperature  and  wind,  respectively.  The  mean-square 
values  of  F  and  G  (see  Eqs.  7  and  8)  at  500  mb  are  denoted 
by  S(F)  and  S(G). 


Cycle 

S(T) 

°C 

S(u) 

kt 

S  (v) 
kt 

S(F)  , 
kt  hr 

S(G)  , 
kt  hr  1 

1 

0.324 

3.83 

4.34 

7.68 

5.75 

2 

0.325 

3.82 

4.02 

6.39 

5.12 

3 

0.323 

3.85 

4.01 

6.20 

5.02 

4 

0.322 

3.88 

4.03 

6.16 

5.07 

30 


TABLE  2.  Ratio  between  observational  and  dynamic  weight  as 
a  function  of  cycle  number.  The  hypothetical  ratios  are  shown 
in  parentheses  and  the  calculations  based  on  values  from 
TABLE  1.  • 


Cycle 

1  (1.5)  ' 

f  (250) 

1 

2.78 

424 

2 

2.18 

327 

3 

2.06 

305 

4 

2.04 

308 

31 


TABLE  3.  The  iterative  solution  to  the  Euler  equations  is 
accomplished  by  the  accelerated  Liebmann  method.  The  number 
of  iterations  required  for  convergence  is  displayed  as  a 
function  of  the  cycle.  These  numbers  represent  the  calcula¬ 
tions  at  700  mb  and  500  mb  for  wind  and  temperature,  respectively. 


3  5  6 

4  4  5 
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APPENDIX:  NATURAL  BOUNDARY  CONDITIONS 


The  right  hand  side  of  (15)  has  categorically  grouped  the 
terms  involving  integrals  over  the  boundary  surfaces.  These 
integrals  originate  from  the  last  set  of  bracketed  terms  in 
(12) .  The  mechanics  of  extracting  the  surface  integrals 
from  the  functional  is  explained  in  Lanczos  (1970).  Before 
writing  the  explicit  form  of  these  integrals,  it  is  convenient 
to  adopt  the  following  notation: 

x',  xM  :  boundaries  of  analysis  space  in  east-west 
direction 

y’,  y"  :  boundaries  of  analysis  space  in  north-south 
direction 

a',  a"  :  boundaries  of  analysis  space  in  vertical 
direction 

Realizing  that  the  global-band  is  cyclic  in  the  longitudinal 
or  east-west  direction,  the  relationship 

4*  (x 1 )  =  <Kx")  (Al) 

holds,  where  <}>  represents  any  of  the  variables,  u,  v  or  T  or 
derivatives  of  these  variables.  The  explicit  form  is: 
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where  [  4>  (x,y,a)  3x,  =  <j)(xM,y,a)  -  4>(x’,y,c) 
notation  of  integral  calculus.  The  terms 
vanish  by  virtue  of  (Al) .  If  6T  =  0  along 
walls,  then  ®,  @  and  (S)  vanish.  When 
upper  and  lower  boundaries,  then  terms  (7) 
These  are  the  Dirichlet  conditions  for  the 
they  have  been  adopted. 


in  conformity  with 
©  ,  @  and  (6) 

the  north  and  south 
6u  =  6v  =  0  on  the 

-  (57)  vanish. 

formulation  and 


300mb,  (S  =  1.25 


ire  1:  Schematic  representation  of  ’’Ascent-Cascade”  process  of  upper-air 
analysis.  The  fields  at  the  surface  and  250  mb  are  control  or  anchor  levels. 
Analysis  begins  at  850  mb  and  proceeds  to  300  mb  (sequencing  from  bottom  to 
top  and  left  to  right  in  the  figure).  After  completion  of  a  cycle  from  850 
to  300  mb,  the  procedure  is  repeated  using  the  adjusted  values  as  input. 


ire  2:  Sea  level  pressure  at  OOZ  on  15  February  1972.  This  analysis  is  the 
result  of  a  two-dimensional  NVA  adjustment  that  has  incorporated  both  wind 
and  pressure  reports.  The  observed  wind  (knots)  and  pressure  (tenths  of  a 
millibar)  rejected  during  analysis  have  been  denoted  by  dotted  barbs  and 
blocked  numbers,  respectively. 
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Figure  3:  Surface  winds  (knots)  obtained  in  conjunction  with  the' pressure  analysis 
is  shown  in  Fig.  1.  This  variational  analysis  served  as  the  lower  boundary 
condition  for  -the  construction  of  the  upper-air  wind  and  temperature  fields. 
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Figure  4:  250  mb  winds  (knots)  which  serve  as  the  upper  boundary  condition  for 

construction  of  the  intermediate  levels.  This  analysis  is  derived  by  the  SCM 
scheme  and  has  not  been  subjected  to  variational  adjustment. 


Figure  5:  The  wind  observations  (knots)  used  by  the  SCM  scheme  to  construct  the 
250  mb  analysis.  The  guess  field  for  the  SCM  above  20N  is  the  12-hr  forecast 
from  the  primitive  equation  model. 
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Figure  6:  ^Vector  difference  between  the  NVA  velocity  and  SCM  velocity  at  700  mb, 

V  -  W .  •  This  vector  represents  the  adjustment  after  two  cycles  of  the  Ascent- 
Cascade  process. 
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Figure  8:  700  mb  wind  analysis  (knots)  by  SCM  scheme,  \V(700mb) . 


Figure  9:  400  mb  wind  analysis  (knots)  by  SCM  scheme,  \V(400mb). 
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Figure  10:  Wind  shear  (knots)  between  400  and  700  mb,  W(400)  -  W(700),  super¬ 

imposed  on  the  temperature  field  at  500  mb.  Observed  temperatures  used  in  the 
analysis  have  been  denoted  on  the  figure. 
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Figure  11:  Wind  shear  (knots)  between  700  mb  and  the  surface  \V(700)  -  W 
superimposed  on  the  temperature  field  at  850  mb.  This  analysis  used  (3 
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